The geometrodynamical origin of equilibrium gravitational configurations 
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The origin of equilibrium gravitational configurations is sought in terms of the stability of their 
trajectories, as described by the curvature of their Lagrangian configuration manifold. We focus 
on the case of spherical systems, which are integrable in the coUisionless (mean field) limit despite 
the apparent persistence of local instability of trajectories even as A*' — !> oo. It is shown that 
when the singularity in the potential is removed, a null scalar curvature is associated with an 
effective, averaged, equation of state describing dynamically relaxed equilibria with marginally stable 
trajectories. The associated configurations are quite similar to those of observed elliptical galaxies 
and simulated cosmological halos. This is the case because a system starting far from equilibrium 
finally settles in a state which is integrable when unperturbed, but where it can most efficiently 
wash out perturbations. We explicitly test this interpretation by means of direct simulations. 
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In a contemporary context, the Neviftonian dynamics of 
a large number of self-gravitating particles describes the 
dynamics of putative dark matter, and to first order the 
formation and evolution of galaxies and clusters. Galax- 
ies that are thought to have formed through largely dissi- 
pationless processes show remarkable regularity in their 
density profiles, and dark matter halos identified in cos- 
mological simulations have similar (and nearly) universal 
profiles [l|, |2| . Yet there is no satisfactory theory predict- 
ing these preferred choices of dynamical equilibria. 

Steady states are characterized by the absence of 
coherent modes, which decay due to differences be- 
tween motions on neighboring trajectories. This effect 
is present in all nonlinear systems (phase mixing [3]), 
but is most efficient if nearby trajectories rapidly diverge 
and spread in phase space. However, attempts at relat- 
ing these trajectory stability properties to the choice of 
A^-body equilibrium — or extracting almost any interest- 
ing information at all from them — are frustrated by a 
puzzling phenomenon. As A^ — ^ cx), particles move in the 
mean field potential, and steady state systems with suffi- 
cient symmetry become near integrable — e.g., spherical 
systems are completely integrable when in equilibrium, 
and phase space distances between their trajectories thus 
generally diverge linearly in time 3]. In contrast, it has 
long been realized that direct integrations of the equa- 
tions of motion invariably reveal exponential divergences 
on a timescale that does not increase with A^ [4, 5]. 

A clue concerning the resolution of this paradox comes 
from studies of softened systems. For when the singu- 
larity in the potential is removed (e.g. 4> '^ \/r — > ~ 
\j\jr'^ -I- const), the divergence timescale does increase 
with A [6, 7]. The purpose of the present study is to 
show that, in this case, physically interesting informa- 
tion, pertaining to the choice of equilibrium, can be ex- 
tracted from the stability properties of trajectories. 



The divergence of nearby trajectories can also be de- 
duced from geometric representations of motion and sta- 
bility in gravitational systems (e.g., Is Mllh . Though 
the geometric description is not unique [1^, one well 
known formulation has the advantage of involving a pos- 
itive definite metric defined directly on the Lagrangian 
configuration manifold — a curved subspace of the 3A^- 
dimensional space spanned by Cartesian particle coordi- 
nates, the phase space being its cotangent bundle |13| . 
This conformally flat metric is expressed in terms of 
the Cartesian coordinates and the system kinetic energy, 

W = E-V,hy ds" = WY.ZN {dx^f (e-g-, 0, l3)- 
The divergence of trajectories is then determined by the 
two dimensional sectional curvatures fcu,n, where u is the 
system velocity vector and n define directions normal to 
it. When the fcu,n are negative, trajectories are expo- 
nentially unstable [15|. Isotropic A^-body systems with 
singular two-body potentials are predicted to tend to- 
ward this state as A^ is increased and direct two-body 
collisions are ignored [SJ, |9| ; and trends similar to those 
obtained from direct calculations are recovered when the 
full gravitational field is taken into account [3|. 

We will be interested in how A^-body systems damp out 
fluctuations, we therefore seek an averaged measure of 
their response to random perturbations. Integrating over 
directions n, we obtain the mean (or Ricci) curvature 
Tu — X]u=i fcn,,,u(s). Furthermore, we will focus on 
isotropic systems. In this case, a further average over 
the velocity directions u can be made, giving the Ricci 
scalar; which, for TV ^ 1, can be expressed as 
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where VW and V^W are the gradient and Laplacian, 
taken with respect to coordinates qi — ^rupXi, and nip 
are the particle masses (which we will assume to be 



equal). For sufSciently large- i V sp herical systems with 
isotropic velocities r„ — >• R/N [16|. 

A correspondingly averaged equation, describing sta- 
bility in the time domain, is given by (e.g., [ll|) 






Q{t)Y^O. 



(2) 



Here Y measures the mean divergence of trajectories and 

Near dynamical equilibrium, W ^ and the dominant 
term on the right hand side will be the curvature term. 
For singular Newtonian potentials, and when direct col- 
lisions are ignored, Eq. ([T]) predicts negative R (since the 
first term vanishes) , and exponential instability is present 
even in the case of spherical equilibria. Ommitng direct 
collisions however implies an incomplete manifold; it is 
not clear if straightforward application of the geometric 
approach remains appropriate in this case [l7|. 

Formally, as A^ — ?> oo, the first term in Eq. ((T]) would 
contribute through a collection of delta functions; but 
if the singularity in the potential is softened, this term 
represent a smoothed out density field integrated over 
the particle distribution (with the softening length act- 
ing as smoothing parameter), and, in dynamical equilib- 
rium, the two terms on the right hand side of Eq. ([1]) be- 
come comparable. Indeed, by multiplying the terms by 
W^/N'^, taking the inverse, and the square root, these 
terms can be seen to be linked to the natural timescale 
of the system (the dynamical, or crossing, time) — the 
first through the relation td ^ l/{GpY/'^ (where the av- 
erage density is related to the Laplacian via the Poisson 
equation), and the second through r^ ^ (W^/mp)^/^a~^, 
where a is the average acceleration affecting a particle. 

In fact, the structure of Eq. ([T]) reflects a pressure- 
gravity balance. To see this, assume i? = 0. Then 
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where pi and a^ refer to the density and acceleration at 
the position of particle i. This equates a density multi- 
plied by the velocities squared (i.e, a pressure term) to 
a gravitational binding term. In the continuum approxi- 
mation, and for a spherical system, this becomes 

/3 /■ m? 

p^r^dr = -MG I -^pdr, (5) 

where p = p{r) is the radial density, m is the mass en- 
closed within radius r, M is the total mass, and the inte- 
grals are evaluated over the volume of the configuration. 
In dynamical equilibrium one can use the Jeans equation 
for an isotropic system (e.g., [3|) to get 
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where a^ = v'^{r) is the one-dimensional velocity disper- 
sion. There is one easily identifiable solution of Eq. ([6]) ; it 
corresponds to p ^ l/r^ and a = constant. This is a sin- 
gular isothermal sphere; though its gravitational force, 
and both terms in Eq. ([S|), diverge as r — >■ 0, and the 
mass diverges as r — > oo, there is significance to this so- 
lution, as isothermal spheres enclosed by a boundary are 
entropy maxima of softened gravitational systems, which 
do tend toward these maxima, provided certain condi- 
tions are met IB, llSJ . 

For well behaved systems, pa^ = as r — >■ and r — > 
cx). Integrating by parts, noting this, and that the kinetic 
energy is related to the average velocity dispersion by 
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\M{a'^)r, and dm = Airr^pdr, one gets 
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pa dm. 
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If well defined pressure and temperature functions can 
be assigned (i.e., assuming local thermodynamic equilib- 
rium), then P = pcr^ and W =^Nk {T)r, and 



P dm = 2 —lJI. I p dm. 



(8) 



Save for the integrations over the mass distribution, this 
form is that of an ideal gas equation of state. Never- 
theless, because of the factor 2, there are no associated 
isothermal states; in contrast, solutions require that the 
density and velocity dispersion distributions correlate. 



APPROACH TO EQUILIBRIUM: AN 
INTERPRETATION 



The results obtained thus far suggest that dynamically 
relaxed equilibrium systems should have small or vanish- 
ing R; but how do they reach such configurations? 

The density and potential of a system started far from 
equilibrium will fluctuate until collective motions are 
damped; this will take place most efficiently when it finds 
a configuration where its response is least coherent; that 
is, when the divergence of neighboring trajectories are 
maximal, and coherent modes are damped through phase 
space mixing of trajectories as the system responds to 
density and potential fluctuations. At the same time, it 
is required that in the absence of fluctuations the trajec- 
tories of a spherical system should not be unstable, so 
as to recover an integrable system in the innite-A^ limit. 
This implies that R cannot be negative at equilibrium. 

Eq. ([2]) suggests that configurations that satisfy these 
requirements have i? — ^ at equilibrium: a system start- 
ing far from equilibrium should then end up in a config- 
uration satisfying this condition. To show this quantita- 
tively, we have used a formulation, developed by Pettini 
and collaborators, which assumes that, for positive Q, 
unstable solutions of Eq. ^ can be estimated by divid- 
ing Q into a mean term fco and rms fluctuations au with 



characteristic timescale r; this gives rise to an effective 
Lyapunov exponent (e.g., Eqs. 79 of [12]) 



x-Ua ^^" 
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^ith A^ = air + ^{^Y + atr^. 

In general, phase space averages are not weh defined for 
gravitational systems (which have a non-compact phase 
space). Nevertheless, every exact {N — >■ oo) coUisionlcss 
equilibrium has a characteristic R and associated fco — 
'^9^N^ - ^^ ^^ therefore possible to quantify the effect of 
fluctuations of a certain magnitude a^ — {Q^) — (Q)^ 
about a given equilibrium, and to compare the magnitude 
of this effect for different collisionless configurations. 

The natural timescale for fiuctuations in a gravita- 
tional system is the dynamical (or crossing) time, so 
T ~ td. And since Q comes in units of inverse time 
squared, we can write fco = a/rf^ and a = b/rl), and ex- 
press the exponentiation time associated with trajectory 
divergence in terms of the natural timescale as 
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This is proportional to ab~^ for a ^ 6 and tends to 
22/3^-2/3 as 6 ^ a. For a given fluctuation level 6, the 
decoherence of trajectories due to local divergence will be 
more efficient, i.e. its exponential time will be smaller, 
for smaller a, and therefore fco and average R. 



TESTING THE INTERPRETATION 

We now examine the above interpretation by means 
of direct A'^-body simulations. Since we are interested 
in estimating the densities and accelerations in the con- 
tinuum limit, we use a technique [20|, which expands 
the density and potential in smooth functional series; 
and since we are interested in strictly spherical config- 
urations, we only make use of the radial expansions in 
this scheme (which is carried up to order 30). Systems 
are sampled using a 100000 particles, started from ho- 
mogeneous spatial initial conditions inside a unit sphere 
of unit mass. This configuration has a natural timescale 
TD = (G/9(0))~^/2 = Y^47r/3, which will be our time unit. 

We have run configurations starting with isotropic ve- 
locities that are either constant or following various de- 
creasing or increasing functions of radius. The results 
shown here correspond to the latter case, for the following 
reason. Collisionless evolution implies that the 'pseudo- 
phase space density' Pp = p/a^ generally decreases along 
particle trajectories fSj; and one of our aims is to explain 
the dynamical origin of cosmological halos, which have (a 
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FIG. 1. Time evolution of normalized curvature Rn (Eg.llip. 
measured in initial dynamical times, for systems starting with 
(from top to bottom) virial ratios, 1, 0.5, 0.25 and 0.125 
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FIG. 2. Equilibrium density (solid line) and phase space 
density (crosses) of system starting with virial ratio 0.125. 
Dashed lines are eye fits withp ~ l/r(l + r- "') and pp ~ r^' 



centrally divergent) pp ~ ^.-i-STS J2l|. From among vari- 
ous functional forms tried, results will be shown for initial 
a ^ r, noting that the trends reported seem generic. 

We vary the initial virial ratio Vir = 2W/\V\ from 1, 
corresponding to near equilibrium, to 0.125, correspond- 
ing to a system started far from equilibrium. Given the 
sampling errors for finite N, R will never tend to zero ex- 
actly in practice; and the absolute values of the densities 
and accelerations (and local dynamical times) strongly 
depend on the central concentrations of the final config- 
urations. A relative measure is thus required for compar- 
ing these. We use the normalized quantity 
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-1, (11) 



which measures the departure from equality in Eq. (4). 

The results are in Fig. 1, which shows that the equilib- 
rium Rn is indeed always of order unity or less, and that 
the further from equilibrium the initial state, the closer 
the evolved system is to a configuration that minimizes 
R. The initial i?„ is zero when Vir = 0.5; systems with 
initial Vir < 0.5 start with R < and then explore state 
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FIG. 3. Evolution of density perturbations, defined as 

J2tM(^ = 0) - P»W)VE^i put = 0), when coordinates 
of final states of systems starting from virial ratio 0.125 (solid 
line) and 1 (dotted line) are perturbed by 10% 



when the singularity in the potential is removed. What 
transpires instead is a description of gravitational equi- 
libria. We propose an interpretation of the approach to 
equilibrium that is effectively a rationalization of the clas- 
sic idea of 'violent relaxation' [23]: systems starting far 
from equilibrium reach such a state once a configura- 
tion where collective oscillations are efficiently damped 
is reached, and such configurations correspond to states 
whose dynamical trajectories are marginally stable. Im- 
posing this condition leads to relations describing relaxed 
dissipationaless equilibrium configurations that are re- 
markably similar to those observed and modeled. These 
systems do not 'ring' and are more robust. 



space until they achieve configurations where (according 
to Eq. [T0|) oscillations are efficiently damped. 

Also, as shown in Fig. 2, the equilibrium density and 
phase space density profiles of this system are remark- 
ably similar to those of cosmological halos (and shallow 
cusp elliptical galaxies [2]). There are some differences. 
Our outer density profile is steeper than the 1/r^ form 
of cosmological halos; however finite configurations can- 
not have 1/r^ outer profiles. (Cosmological halos are not 
isolated equilibrium structures as their outer regions are 
subject to mass infall.) The functional form of our empir- 
ical spatial density fit is also slightly different from stan- 
dard 'NFW' (where pnfw ~ r(i+r)^) )' ^^'^ there seems 
to be some flattening toward the center of the l/r cusp. 
Yet, similar departures from ideal NFW often also apply 
to cosmological halos [22]. 

Finally, we explicitly verify that configurations with 
smaller i?„ do efficiently damp out fluctuations. This is 
already apparent in Fig. 1, where uctuations are quickly 
damped for such systems. (This may be intuitively ex- 
pected, for these systems lack homogeneous harmonic 
cores with nearly constant orbital frequencies and co- 
herent response.) Further verification is obtained by 
perturbing the final states and measuring variations in 
the density distributions as systems, so perturbed, are 
evolved. The results of one such experiment, where the 
spatial spatial coordinates were decreased by 10%, are 
shown in Fig. 3. There, the system that was perturbed 
from a state with smaller i?„ is seen to quickly settle 
back to an equilibrium that is closer to the unperturbed 
one, and to gradually close in further on it (as the outer 
regions, with large crossing times, evolve). 



CONCLUSION 

Contradictions between predictions pertaining to the 
dynamical stability of large- A^ body system trajectories 
and those implied by the collisionless limit disappear 
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